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ABSTRACT 



The general differential equations describing unsteady-state heat 
transfer with a fluid flowing through a porous medium are derived. These 
equations represent a physical model for heat transfer in thermal oil- 
recovery process, packed-bed chemical reactors, and heat regenerators. 
Fluid-solid convective heat transfer and longitudinal conduction in both 
the fluid and solid phases are considered. Laplace transformation and 
numerical inversion are used to solve the system of partial differential 
equations. A digital computer program obtains the numerical results 
which are compared to those of Green and Perry using finite difference 
technique, and to experimental data of Preston. Also presented are 
analytical solutions for the cases where the longitudinal conduction is 
neglected and the convective heat transfer coefficient is assumed to be 
infinite. These solutions are programmed and results are compared to those 
from the general case. The effect of different heat transfer mechanisms 
on temperature profiles at low fluid velocities is studied. The results 
show that this numerical method gives accurate results with relatively 
short computational time. 
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NOMENCLATURE 



Symbol 




Quantity 


Unit 


a 


= 


Surface area of solid particle per unit 
of bulk volume 


ft2/ft3 


A 


= 


Total heat transfer area 






= 


Fluid cross sectional flow area 


II 




= 


Matrix cross sectional flow area 


II 




= 


Fluid phase specific heat 


Btu/lb.°F 




= 


Solid phase specific heat 


II 


Cl 


= 


0 


II 




= 


1 

+ 


II 




= 


Average particle diameter 


ft 


F 


= 


Number of time units per hour 


UT/hr 


h 


= 


Heat transfer coefficient 


Btu/hr ft°F 






Effective thermal conductivity of porous 
medium, assuming solid and fluid temperature 
are equal 


Btu/hr ft°F/ft 


k 

s 


= 


Effective thermal conductivity of solid phase 


II 


K 


= 


k, (1-0) 


Btu/hr ft^oF/ft 




= 


Effective thermal conductivity of fluid 
phase 


II 


k ' 


= 


k^0 


1 1 


'^fc 


= 


Molecular thermal conductivity of fluid 


II 


^fm 




Effective conductivity of fluid phase due 
only to fluid mixing or dispersion in 
porous medium 


II 




= 


Static effective thermal conductivity of 
porous medium 


II 


L 


= 


Length of packed bed 


II 


s 


= 


Laplace transform variable 


dimensionless 



Symbol 




Quantity 


Unit 


t 


= 


Dimensionless time parameter, ha0 


dimensionless 


Tf 


= 


Fluid temperature 


OF 


Ti 


= 


Injected fluid temperature 


Of 


T 

s 


= 


Solid temperature 


OF 


u 


= 


Solid temperature fraction, I's^'^i 


dimensionless 


V 


= 


Fluid temperature fraction, Tg/T^ 


II 


Vf 


= 


Fluid interstitial velocity 


ft/hr 


UT 


= 


Time unit 


fraction of h: 


X 


= 


Distance from point of fluid injection 


ft 


X 


= 


Dimensionless distance, x 

L X 


dimensionless 


Y 




Dimensionless distance, ^ha y 


tt 






Fluid mass flow rate 


Ib^/hr 


Ws 


= 


Mass of solid matrix 


B 

r— 1 


N 

re 


= 


Modified Reynolds number, 


dimensionless 




= 


Peclet number ^1* 

0< 


II 


NTU 


= 


Number of heat transfer units hA 

cbfCf 


It 



GREEK LETTERS 



c< 

P 

P' 

r 

A 

A' 

d 



= Thermal diffusivity, 






= Ratio of thermal dif fusivities , 



0(: 






ft^/hr 



dimensionless 



= Dimensionless time 
= Time 



■ fef 








II 




II 




II 


■ , As 


It 




II 

hr 



vii 



Symbol 

P 

0 



Quantity 

Viscosi ty 
Density 

Porosity of porous medium 

Ratio of heat capacities per unit length 



Unit 

Ibj^^/ft hr 
dimensionless 



NOTE: An occasional term may appear in the body of the text that does 

not appear in this list. Such terms are used only once and are defined 
as they appear. 
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1 . 



Introduction . 



Thermal recovery operations are rapidly growing in importance 
throughout the oil producing industry. Large volumes of oil previously 
considered uneconomical to recover are being produced by thermal processes. 
The intense interest in the application of the thermal energy to oil 
reservoirs as a means of increasing the percentage of oil recovery has 
stimulated the research on the problem of heat transfer in porous media. 
Thermal recovery has seemed most applicable to reservoirs that contain 
very viscous oil at reservoir conditions. This is due primarily to two 
factors: the low recovery from viscous oil reservoirs by primary pro- 

duction or conventional secondary recovery methods and the significant 
reduction in viscosity that takes place when viscous oil is heated. 

In these thermal methods, heat is injected or generated in the 
reservoir. The heated oil has its viscosity decreased thus making the 
removal from the reservoir easier. Thermal energy may be applied to a 
reservoir in different ways. The simplest processes are steam injection 
and hot water injection. In more complicated processes, the crude oil 
is burned at one end of the reservoir, forming a combustion zone which 
moves toward the other end. The product of combustion is a mixture of 
oil and condensed water, resulting from thermal cracking. No matter which 
method is used, the effect of heat on the production of oil and water 
should be known. Thus, a knowledge of various heat transfer mechanisms 
with their individual effects is required and also the temperature 
history at each point of the reservoir and the temperature distribution 
throughout the reservoir should be known. 

Previous studies on heat transfer in porous media may be classified 
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into two groups. The first group considered that the main heat transfer 
mechanism involved in this problem is convection from the fluid surface 
to the solid surface. Thus the longitudinal conduction in both the fluid 
and solid phases are neglected. This case was intensely studied by 
Anzelius [l ] , Hausen [lO] , Nusselt |^22] , Schumann [ 29 J , and 
others. The second group assumed that the film resistance is negligible 
and that the heat is transferred solely by longitudinal conduction. This 
attack on the problem was made by Jenkins and Aronofsky [ 14 ] . Preston 
[ 24 ] used their solution to compare with the results from his experi- 
mental work. Authors on the problem of heat regenerators considered 
the matrix of the heat exchanger as a porous medium through which a 
gas is pumped. In this case the stored energy and the longitudinal 
conduction in the fluid were neglected. Green and Perry [?J have 
investigated the general case where both conduction in the direction of 
flow and convection from fluid to solid were considered in the mechan- 
ism of heat transfer. They used finite difference techniques to solve 
the general set of partial differential equations. This forward dif- 
ference approach has its disadvantages because of the small time and 
space increments necessary. 

It was the purpose of this thesis to use Laplace transforms to 
solve the differential equations derived from the heat balance for the 
general case. A FORTRAN program was set up for use with a GDC 1604 
digital computer. By using the parametric values of Green and Perry |^73 
and of Preston [^24]) , numerical solutions were obtained and checked 
against their solutions. Also, analytical solution to the differential 
equations of Schumann and Jenkins-Aronof sky were derived by using 
Laplace transforms. Two programs were set up and numerical solutions 
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were compared with the solutions to the general set of equations. The 
purpose of this comparison was to determine the relative importance of 
the different heat transfer mechanisms and how these mechanisms are 
affected by changes in the significant parameters. 

The general case studied in Section 3 is limited to a model of 
infinite length. The mathematical derivations for the case of heat 
regenerators and of packed beds of finite length are presented in 
Appendix III. An outline is presented here of additional work which 
would be required to produce numerical results for this case. 
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2 . Literature Survey 

The theoretical solution to the problem of transient heating of 
porous media should provide: 

a. The temperature history at a point in the porous medium as 
a function of time. 

b. The temperature distribution throughout the length of the 
medium at a given time. 

The mathematical and physical model of an oil reservoir is similar to 
that of a packed-bed or of a heat regenerator. Many theoretical studies 
have been made in these areas. The reservoir can be considered as a 
semi-infinite porous body through which the fluid is flowing. The fol- 
lowing assumptions are usually made: 

a. The initial solid and fluid temperature are equal throughout 
the length of the body. 

b. The fluid is injected at one end. At time zero, its tempera- 
ture is suddenly changed to a higher value and kept constant at this 
end , 

c. The rate of fluid flow is constant. 

d. The physical properties of fluid and solid are independent of 
temperature • 

e. No temperature gradient exists in the direction perpendicular 
to the flow direction, i#e., the conductivity of the solid is infinite 
in that direction. 

The basic mechanisms of heat transfer in a porous medium through 
which the fluid is flowing are: 

(1) Storage of heat in an element of fluid. 

(2) Conduction of heat through the solid and the fluid phases. 
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(3) Convection between the solid and fluid phases. 

(4) radiation. 

Radiation may play a significant role in the energy transfer encountered 
in the problem of transpiration of fluid in chemical reactors, heat 
shields and solar heat collectors. This mechanism is assumed negligible 
in an idealized model of a thermal oil recovery process, packed-bed 
chemical reactors and heat regenerators. The differential equations 
applied to the general case where both conduction and convection are 
considered can be derived from heat balance as presented in the next 
section. The original equations are: 

For the fluid phase : 

Ti- - (1) 

For the solid phase : 
where : 

Subscript f refers to fluid phase . 

s refers to solid phase , 

T = temperature above a base temperature which is the initial bed 
temperature, ^ 

2 o 

a = surface area of solid particle per unit of bulk volume, ft /ft*^ 
c = specific heat, Btu/lb.°F 

2 o 

h = heat transfer coefficient, Btu/hr.ft . F 
k = pseudo-thermal conductivity, Btu/hr . f t^ . °F/f t 
X = distance from point of fluid injection, ft 
= linear velocity of fluid, ft/hr 
0 = bed porosity, dimensionless 

9 = time, hours 

= density, Ib/ft^ 



P 
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The different mechanisms of heat transfer involved in the general heat 
balance equations were discussed in detail by Hadidi []9] . 



Since an analytical solution to this set of differential equations 
is obviously difficult, all previous studies were confined to special 
cases, where either conduction or convection is neglected. An outline 
of this literature might be helpful to the reader. 

Case 1 : k = 0, 0<^Ka<<t>o'' 

By assuming that conduction in both phases is negligible, one can reduce 
the equations (1) and (2) to; 



alternative method of solution to the same equations in a thermal 
analysis of the counterflow regenerative heat exchanger. Nahavandi 



to present a solution to the rotary heat exchanger problem. Lambertson 



of a thermal regenerator by using finite differences to solve this 
problem. 

In Appendix I the writer derived the solution to this special case 
by simply using Laplace transforms. The same dimensionless parameters 
in the general case were used and the solution was then programmed to 
provide numerical data which were compared with the solution to the 
general problem. An alternate attack to the problem was made by 
Creswick [5] • In his analysis, he neglected the term p^c^0 




(4) 



(3) 



Hausen 




etc. Several techniques have been developed to solve the 



system of equations (3) and (4). C. E. Iliffe 




presented an 




used Laplace transform and power series expansion 




and, recently , A. J. Willmott presented a digital computer simulation 




6 



in equation (3) describing the heat gained by an element of the moving 
fluid, but he considered important the effect of longitudinal conduction 
in the solid by adding the term k^(l-0) ^ to equation (4). These 

t> 

two equations were solved by finite difference techniques. Bahnke 
used Creswick*s equations and finite differences to solve for the con- 
duction effect on effectiveness of the rotary regenerator. Recently, 
Moreland [^20] applied Laplace transform and Gaussian quadrature for 
numerical inversion to get the solution to the “single blow” problem, 
using the same set of equations. 

Case 2 : 0<k<oo 'ha=oo. 

In this case, the fluid-solid boundary resistance was negligible, i.e., 
ha is infinite or T£ = T^, But the porous body is considered as a 
homogeneous unit with longitudinal conduction in the direction of flow. 

By substituting for the term ha (T£ - Tg) in equation (1) its 
value derived from equation (2) and letting T£=Tg=T, we have a combined 
equation: 

where r: effective thermal con- 

ductivity of the porous medium. This approach to the problem of heating 
of porous media was offered by Jenkins and Aronofsky [^14j . The writer *s 
solution to equation (5) was derived by using Laplace transforms and 
is presented in Appendix II. A program was set up to provide numerical 
results which were checked against the solutions to the general case. 

Jenkins and Aronofsky, after investigating the results and checking 
them against published data, mentioned that by selecting a value of k^ 
which gives the best agreement between experimental temperature profile 
and the analytical solution to equation (5) , one can determine the com- 
bined dynamic thermal conductivity of the porous system. 
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Preston [24] in his experimental work, measured the static thermal 



conductivity of the porous system under no-flow conditions and the 
dynamic thermal conductivity under flow condition. He concluded that 
at velocities less than 0.05 ft/hr, (i.e., at velocities characteristic 
of flow in petroleum reservoirs), the effective thermal conductivity 
would equal the static thermal conductivity. For greater velocities, 
he stated that the effective thermal conductivity under flow condition 
increased with velocity. He concluded that at low rates of flow, the 
main mechanism of heat transfer in porous media could be assumed to be 
longitudinal conduction alone, i.e., the convection heat transfer 
could be neglected. 

Case 3 : Both k and ha are finite 

This is the most general case where both conduction and convection are 
assumed to be important. The general set of differential equations 
given by equations (1) and (2) is too complicated for an analytical 
solution. In a recent paper. Green and Perry [^7] reported the results 
obtained from a numerical analysis of the problem. They reduced these 
equations to difference equations of the forward difference type and 
solved for several values of parameters on an IBM 650 digital computer. 
Their solutions were checked against the results of Preston [^24] with 
close agreement. The writer’s approach to solve the system of dif- 
ferential equations (1) and (2) is presented in the next section. A 
computer program was set up for use with a CDC 1604 computer. Numerical 
solutions were compared with the results of Green and Perry [^ 7 ]] and of 
Preston [24] . 
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3. 



Mathematical Analysis 



The differential equations (1) and (2) can be derived from heat 
balance as follows; 

a. For the fluid phase : 

Heat stored in an element of fluid = C 0 



Convection by moving fluid 

Conduction in the fluid 

Heat transfer to the fluid element 



h'i ^ 



0 



by convection = 

The heat balance yields: 

This is the same as equation (1). 



b . For the solid phase : 

Heat gained by an element of solid = 

Heat transferred to the solid 

element by convection = ha ( Tf - ”^s) 

Heat transferred by conduction 

from the solid element = U. Ci_0') i- 

The heat balance gives: 



( 6 ) 



+ HaCT^-T,) (7) 

This is the same as equation (2) 

Let us define two new dimensionless variables: 



= dimensionless 


distance 


_ 










vkfo; 


= dimensionless 


time 





f— t-v.e 








VkfOj + 



Substituting these new variables into 

Hi =-HL ^ f f 



(6) , we get: 
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We introduce the dimensionless parameter: 

A = where «, = 

Equation (8) becomes: 



7>lf 


:>Tt 

■ ■ 


- -MTf-T.) 


(9) 


By the same substitution 


, equation (7) becomes: 




ST* 

2 T 


D ^ 




(10) 


where 


C<|. = 


ks 

P. c. 














i<; = 






Let U, = 


and -yr - -ii- 

■’"t Tc 




where is the injected 


fluid temperature 




II 


— , and -8 


-hi 

K's 




the system 


of equations 


is then: 




7>^ 


= (Afi) ^ 


+ ()ip)(v-a) 


(11) 


li X 




- _ ACv'-^) 


(12) 



Boundary and initial conditions 

initial fluid and matrix temperatures are uniform and 
equal. The base scale temperature can be chosen so that these 
temperatures can be taken as zero for convenience: 

u.(^,o) = v(<i.,o) = o • 

(2) - At y = o and T = o"^, the injected fluid temperature is 

suddenly changed to a different higher value and held constant there- 
after. The input temperature is thus a step function: 

lr(o,^) = i 
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(3) ^ - At y = 0 and T = the solid temperature is assumed to 
instantaneously rise to the value of the step input temperature of 
the fluid: 



^(o,r) = i 

(4) - As y approaches infinity, for all T , the fluid and solid 

temperatures decrease to their initial value; 

Kju{oc,,X) =. \r(oo,‘t) = o 

By calculating tc from the equation (12) , and differentiating it with 
respect to 'C and ^ , we get the following equations: 



IL 


1 

X 


1>X 


A + 

a ^ 


b w" .Xu- 

J 




1) tt 

D r 


1 

~ X 


^ z 




1 Z\ 
X 


() V 


dV 


1 

~ X 




3 


1 

X 





(13) 

(14) 

(15) 



Substituting the terms u, , and into equation (11) and 

It 

rearranging the terms yields the following differential equation in tr : 



(sy) Ir 4- + (&r + i):^ + _L^ = 0 

Let IT be the Laplace transform of ir and s the transformed variable: 



(16) 






The following formulas for Laplace transformation are used: 




XT 







- sxr - ir(i^,o) 



(from the initial condition) 
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If (^,0) 




Substituting for the partial derivatives in the equation (16) their 
Laplace transforms yields the following subsidiary equation: 



A X ^ _ 

I d > ct 






f + sj 



ir 



(17) 



This is an ordinary differential equation in V 
auxiliary equation is: 



The corresponding 



a) r** - p r^_ + +^A(v+l) r^+ (|i'X'+^)r4 



[x 



0 



(18) 



The general solution in the Laplace s plane for the fluid temperature 
is then: 






= G^(s) e’’’^ 4 - Gj(s)e*"^-h C^(s)e.^^+ G^(s)e 



’'•.a- 



(19) 
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where To. / 1^5 / are the roots of the equation (18). 

Since some coefficients of equation (18) are functions of s, which is 
a complex number, the roots ^re expected to be 

complex numbers. 



The boundary conditions (2), (3) and (4) are transformed and then 
used to determine the constants Ga • 



BC 


(2) : 


ir 


(0,0 = i 


BC 


(3) : 


U/ 


(o,s) = ^ 


BC 


(4) : 


d (00 ^ s) 


= v(oo,s)= 0 


Applying 


the BC 


(4), it is 


observed that may only be zero if 



the exponents in equation (19) are negative, i.e., if the real parts 
of the roots are negative. Since the parameters in the coefficients 
are unknown, the number of roots with negative real parts cannot be 
predicted by using Routh’s criterion. Many computer test runs were 
made to investigate the behavior of the roots by using a wide range of 
parameters. Results from these tests showed that only two roots have 
negative real parts. The following derivations were based on that 
remark. If and are assumed to be the roots with positive real 

parts, then and in equation (19) must be zero and \r is reduced to: 

V = Cj^(s)e'"’^+ Cj,(s) (19a) 

Applying the BC (2) to the above equation, one obtains: 

+ C^Cs) = ^ (20) 

Taking the partial derivatives of given by equation (19a) 

and evaluating them at y = 0 gives: 

(0, s) = Gj (s) + G^ (s) 

^ (o,s) = G^(s) + G^(s) 



% 
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We calculate the transform of equation (13): 



UU 






( 21 ) 



Applying the BC (3) to the above equation and using the values of 

— 

and ^ (q above . we have: 

h'--' ^ 

a(o,s)= ^ C„W * 4.^r.C„(s)^. 



or 



/TV = 1 

£<C,Cs) _ i Jr,G, 



1 



Let : 



Ur 1 T 1 r 1 

then equation (22) may be written as: 



n = dL 



= t 

, n. = 1,2 



±_ 

>\ 



(22) 



(23) 



The functions and may be found from the matrix equation; 

_±_ 

S 



1 

Z, 



1 

z. 



Gx 



A 



It follows that: 



Gi = 



Cx = 



Za 

S 



± 

Jl 



2.- 



i 

A 



Z, 

s 



Zx- 



Substituting and into equation (19) results in: 






2 1 



fr-vK* fr-t) 






The value of tv(^,s) can be given by equation (21) : 



(24) 






= _ ^ 
"b u.^ 



— 1 
A 2)^ ^ A ^ 



= ^£-.C„e'-V 

^=1 'ITSI 
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To determine the values of the fluid temperature and solid temperature 
at any given time and at any point of the bed, the inverse transform 
of equation (24) and (25) should be found. This requires the evaluation 
of the inverse integral: 

c - Loo 

Where C is a real constant that exceeds the real part of each of the 
singular points of F(s) which is v(^^s) and in this case. 

If the analytical form of the function were known and 

its poles and branch points could be located without difficulty, the 
inversion integral in (26) could be evaluated by a suitable deformation 
of the path of integration and the use of Cauchy's theorem on the 

residues. But in our case , U and tC (^y^) aire functions of the 

r 

roots of the quartic equation (19) . These roots could be analytically 
calculated but their analytical forms would be so complicated that the 
evaluation of the inverse integral is hopeless. Thus, given numerical 
values of the parameters, and could only be evaluated 

as functions of s. Then some approximate numerical inversion scheme 
must be used. 

The need for inverting Laplace transforms has been experienced in 
many fields, and approximate inversion methods have been developed in 
connection with several subjects. Thomas L. Cost [4] in applying 
numerical Laplace transform inversion to viscoelastic stress analysis, 
presented a unified treatment of the most promising approximate in- 
version methods in his paper. Among these, the orthogonal polynomial 
inversion methods of Papoulis (^23 J| and of Lanczos jj.8^ are mathematic- 
ally well founded. Legendre and Laguerre orthogonal polynomials were 

used. Recently, Moreland [~20^ in his thesis on the "single blow" 

15 



problem, used the technique developed by H. Hurwitz for numerical 

quadrature of Fourier transform integrals and adapted by L. S. 
Schmittroth for the inversion of Laplace transforms* But the 

most sophisticated investigation on the numerical inversion method was 
made by Salzer . In his early paper in 1955, he derived the 

properties of a certain set of orthogonal polynomials > that 



play a role in inversion integrals 



ioo 




, similar to 



those of the Laguerre polynimials which %ve used to evaluate the direct 

oo 

Laplace transform integrals . A short table of 



weights and zeros was also furnished by the author. In his later paper 
, Salzer presented a table of weights which may be used in con- 
junction with values of F(s) evaluated at integral values of s, using 
the formula: 



C -*• 0 oO 



iiti 




F(s) ds 



(m) 



(t) F(k) 



k = i 



C - iob 

This method is suitable for hand computers since both the weights and 

zeros are real numbers. In another paper on Laplace transforms [27j , 

he presented an extensive table of complex zeros and Christoffel numbers 

up to order n = 16 and with 15 significant figures, for use with his 

first method. Since this is the most promising method, it has been 

chosen to invert the ir(y.,s) and given by equation (24) and 

(25). An outline of this method is presented here. 

Salzer stated that if F(s) is really the Laplace transform of a 

function f(t), it must behave like a polynomial in the variable 1 

s 

without a constant term along the line C-too^ Ci-loo , Then one may 
find f(t) numerically using new quadrature formulas similar to those 
employing the zeros of Laguerre polynomials in the direct L.T.or the 



16 



zeros of Legendre and Chebyscheff polynomials in the methods of Lanczos 
and Papoulis. Suitable choice of yields an n-point quadrature 



formula that is exact when is any arbitrary polynomial of the Irv 

degree in dc = ^ without a constant term. Thus: 
s 

tv 

j P- (i-) P« ij-) 



(26) 



In the above formula, = — are the zeros of the orthogonal 

polynomial b^(x) ^ ~ ^k) ^arived from the generalized Bessel 

k = i ^ 

polynomial defined by H. L. Krall and 0. Frink |[l6 J as: 

dl (^] 

ds"'\ S'^Jf 



kaO 



PtvCy) proved to be 1) e S ^ 

The orthogonal property is: 

C-f loo 



— . r P'rt (4-) 



= 0 



k = 1 ,2, . . . ,n. 



{ix) ^ ** UoO 

in formula (26) are Christoffel numbers. 



There is no loss of generality if equation (26) is written as: 

k = l 



C,- Loo 



where u = st • 



so that F^ujis still a polynomial in 1^ , without a constant term, if t 



is specified numerically. The roots 






and the Christoffel numbers 



a(h) 

are all complex, except when n is odd. They were provided in 
Salzer's paper [^27 J . 

Recently, N. Skoblia of the Academy of Sciences of USSR Moscow, 
has published a booklet |^3(^ presenting tables for the numerical inversion 
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of Laplace transforms. This method is more general than Salzer*s 
method. It enables one to evaluate: 

c 

c- coo 

where C>0 and C lies to the right of all singularities of F(s) and 
m = 0.1(0, 1)3.0. The case m = 1 has been treated by H. Salzer in his 
papers referred to above. The difference between these two methods is 
that Salzer *s quadrature formula is exact if F(s) is a polynomial in 
of degree 2n such that F(oo ) = 0. In Skoblia*s method, the quad 
rature formula is exact if F(s) is of degree (2n-l) , but F( oo ) need 
not vanish. Thus the Christoffel numbers in Salzer *s table differ from 
those of Skoblia but the zeros are the same. Since these tables are 
not yet available at the USNPGS Library, a comparison of the results 
from two methods has not been possible. 
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4 . Computer Programming 



Preliminary test programs were set up to investigate the behavior 
of the roots of the auxiliary equation (18) . Three Library subroutines 
solving the polynomials with complex coefficients were used. The sub- 
routine R00TS2 using MULLER *s method proved to be unsuitable for this 
unusual equation. The roots did not converge after 25 iterations. 
Increasing the maximum number of iterations from 25 to 50 produced con- 
vergence but the computational time was too long as compared with the 
other subroutines. The subroutine COMSUB using Newton Raphson method 
was the fastest but for some range of data, it failed. Finally, the 
subroutine POLYRT using LEHMER's and NEWTON *s method was tested. This 
subroutine gave satisfactory results although it was still slow. A 
combined program using both POLYRT and COMSUB was set up. Since the 
zeros of Salzer*s polynomial do not largely vary from one to another, 
for each set of coefficients, the first run used POLYRT; the roots 
provided by that run are used as guessed roots in subroutine COMSUB. 

Each time the COMSUB fails the POLYRT is used again. 

The subroutine VUBARl corresponds to the general case. It cal- 
culates the roots of equation (18), selects the roots with negative real 
parts, computes and C 2 and complex quantities it ; XT , then sends 
them back to the main program TEMFLM with each set of zeros and Christ- 
offel numbers corresponding to an order m. This enables the main pro- 
gram to provide a set of values of \T and IL . It was found that in- 
creasing the order of polynomials in (-!-) did not necessarily increase 
the accuracy of IT or it . Plotting the values of \T vs the order m 
showed that IT does not converge as m increases. On the contrary, the 
values oscillate at random. Attempting to choose an optimum m also 
failed. But all values of XT agree to at least 3 significant figures. 
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Test runs for the limiting cases showed that the numerical results 
obtained by numerical inversion agree with the analytical solutions 
(provided by the programs Schumann and Jenkins) also up to 3 decimals. 
Finally, in order to shorten the computational time, it was decided 
that only the zeros and weights of order from 11 to 16 were used and 
the average of six values of XT and u. was taken. For each curve of 
XT vs distance or time, only a limited number of points were calculated 
from numerical inversion. The intermediate points were interpolated by 
the subroutine AITKENF . 

The analytical solutions for special cases in Schumann and Jenkins- 
Aronofsky problems were derived by the writer and presented in Appen- 
dicies I and II. Their numerical solutions were provided by programs 
SCHUMANN and JENKINS. Results from these programs were compared with 
the results of the general case to show the relative importance of 
various heat transfer mechanisms. 

The mathematical derivations applied to the case of heat regener- 
ators and of packed-bed of finite length are presented in Appendix III 
where the usual dimensionless parameters in heat regenerator problems 
were used. The problem turned out to be more complicated because all 
the roots of the characteristic equation must be used. Additional 
boundary conditions were used in order to be able to calculate all the 
four constants G(s). A subroutine may be written for this case and 
numerical results may be compared with Creswick*s results [s] . 
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Fig. 1 - Fluid-solid temperature differences; 

effect of dimensionless parameter X 
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FLUID TEMPERATURE FRACTION 



Figure 2. Comparison of generalized numerical solution to 



simplified analytical solution; effect of dimensionless 
parameter A . 
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FLUID TEMPERATURE FRACTION 



Figure 3. Comparison of generalized numerical solutions to simplified 
analytical solutions; effect of solid phase thermal 
conductivity. 
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FLUID TEMPERATURE FRACTION 



Figure A. Fluid temperature profiles; effect of 
dimensionless parameter A , 
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FLUID TEMPERATURE FRACTION 



Figure 5. Fluid temperature profiles; effect of solid 
thermal conductivity . 
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Figure 6 



Comparison of numerical temperature profile to PRESTON'S 
experimental data . 
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Figure 7. Comparison of numerical temperature profile to 
PRESTON'S experimental data. 
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Discussion of Results 



Figure 1 shows temperature profiles for different values of 
dimensionless parameter A • Different values of ^ are used to prevent 
the curves from falling on top of each other. These curves agree very 
well with the results of Green and Perry's calculated by finite dif- 
ference methods. It is seen that as A increases, the temperature lag 
between the two phases decreases. This point can be easily explained 
by the fact that the dimensionless parameter A is proportional to the 
heat transfer coefficient (ha)^ and to the fluid thermal conductivity, 
and is inversely proportional to the fluid velocity. Thus, the temper- 
ature lag between the two phases decreases for large values of ha or 
kf and at low fluid velocities .Thus one can conclude that at very low 
flow rate, as in the case of oil reservoirs, the approximation T^ = Tg 
is reasonable; on the contrary, it is not applicable to the case of heat 
exchangers where the gas flows at high velocity. For very large values 
of ha, one can assume that the fluid and solid temperature are equal. 

In this case, the general partial differential equations are reduced to 
one equation in T. 

The effect of the heat transfer coefficient on the temperature lag 
can be easily studied in Figure 2 where results from numerical solutions 
to the general differential equations are compared to results from 
simplified analytical solutions using the equation of Jenkins and 
Aronofsky and of Schumann. The writer's results also agree with the 
results of Green and Perry. The graph shows that for values of A equal 
or less than .114, the curve for ha, ks and kf finite approaches the 
curve obtained by neglecting the thermal conductivities kg and k^ . For 
values of A larger than .342, the numerical solutions become closer 

to the solutions based on the assumption that only the conduction is the 
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important heat transfer mechanism. The curves also show that in the 
intermediate range of A , both conduction and convection are important 
and should be considered. 

Figure 3 shows the effect of thermal conductivity of the solid 
phase. For values of ^ equal or less than .883, both the simplified 
analytical solutions are acceptable. As the ratio increases, the 

assumption of T£ = Tg is increasingly better. 

Figure 4 shows the effect of dimensionless parameter A on the 
time- temperature history. Decreasing ha or increasing fluid velocity 
makes the temperature at a given point more responsive. The same results 
can be obtained by increasing the solid thermal conductivity. 

In order to compare numerical solution of the general case to 
experimental results, the data of Preston was used. The parameters 
needed for calculation are the thermal conductivities kg and k^, the 
densities ^^d , the specific heats Cg and Cf, the porosity 0 

the heat transfer coefficient ha. It should be pointed out that kg is 
not a true but pseudo-thermal conductivity characterizing the rate of 
apparent solid phase conduction and that k^ is the sum of two terms: 




where kf^ is the fluid molecular conductivity and k£^^ is the term 
characterizing the effect of the eddy-dispersion. This dispersion 
effect is due to the irregularities of fluid flow in packed beds 
causing a convec tive*mixing process. The values used for the parameters 
were furnished by Preston’s data ^24^ , except kg and k£ . In his 
experimental work, Preston measured the thermal conductivity k^, of 
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the porous medium under no flow condition which he called static 
thermal conductivity, and the effective or dynamic thermal conducti- 
vity of the porous system under flow condition. Using this concept 
of static thermal conductivity. Green and Perry [?] assumed that this 
conductivity is the sum of two terms independent of the fluid velocity: 

Thus kg could be calculated from this relation, knowing the values of 
kfc, k° and 0 

The mixing term of the fluid conductivity could be calculated by 
two methods: 



(1) By assuming that the heat transfer Peelet number Npg is defined 



as 









and equal to the mass transfer Np^ , one can use the 



plot of Npg vs available in mass transfer experimental work in 
porous media [6"| to get 

(2) By using the correlating equation derived by Green and Perry 
and Babcock from experimental data [sj : 



Jliin. - o.il5( ] 

kfc V D / 

where D is the molecular diffusivity for fluid phase heat 

transfer. 

Since the reference [6j was not available in time, the correlating 
equation in method (2) was used to calculate The result is not 

reliable, for Green and Perry state that this equation should be applied 
only to the values of greater than 0.03. 

Temperature history was plotted for two systems of packed bed. ,Fig. 6 
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shows that numerical results agree with experimental data. In Figure 7, 
the solid temperature curves approaches the experimental data while the 
fluid temperature curve is not too close. This might be due to the fact 
that the value of this case is beyond the range for which the 

application of the correlating equation is valid. 
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7, 



Conclusions 



The following conclusions may bo drawn from the results discussed 
in the preceding section; 

(1) Approximating the fluid and solid temperatures by the same 
value is reasonable only in the range of very low fluid velocities. 



ties less than 0.05 ft/hr, the effective thermal conductivity under flow 
conditions is equal to the thermal conductivity of the system measured 
without fluid flowing. 

(2) The approximation of T£ = Tg is still applicable to the porous 
systems which have large heat transfer coefficients. 

(3) The fluid velocity considerably affects the temperature profiles. 

(4) At high rate of fluid flow, the heat transfer coefficient plays 
a predominant role (it increases with velocity) . 

(5) Salzer^s method of numerical inversion of Laplace transforms 
may be very helpful for the solution of a system of partial differential 
equations with constant coefficients. It was shown that this method is 
much faster than other numerical inversion methods using Legendre, 

Laguerre and Chebyscheff polynomials; it has also the advantage over the 
finite difference methods which require small increments in space and in 
time for acceptable accuracy. The only problem encountered in this 
method was that the convergence of the integration could not be obtained 
as the order of polynomials was increased. But the results were con- 
sidered satisfactory since values of the inverse transform agree to three 
significant figures for all orders from 4 to 16. Numerical inversion 
methods of Laplace transforms are still in development and promise to be 
the main alternative to the finite difference method. 



This confirms the conclusion of Preston 
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8 . 



Recommendations for Future Studies 



(1) Green and Perry suggested that an effective thermal conductivity 
of a porous system can be obtained by selecting a value of which will 
give the best agreement between the numerical temperature profile and 

the analytical solution to the equation considering Tf = T^. This may 
be done by comparing the maximum slope of the two curves. A prediction 
of the behavior of the slope might be helpful; it could be made by 
examining the formula giving the slope of equation ( 5 ) based on the 
assumption Tf = Tg . This formula derived by Preston [^24]jis: 

dv ^ X e ^ ® 

d 0 aoyii KO' 

where K 5 : ^ 

and kg is the effective thermal conductivity of the system. It is 
shown from the slope formula that the slope decreases as ke increases. 
Results from this recommended investigation may be compared to experi- 
mental data of Preston and then may verify the validity of the suggestion 
of Green and Perry. 

(2) Skoblia*s new method of numerical inversion of Laplace trans- 
forms [ 30 J may be used to solve this problem and make comparisons of 
the two methods. It is expected that Skoblia*s method will give more 
accurate results. 

(3) A subroutine may be written for the case of heat regenerators 
and packed beds of finite length, using the mathematical derivations in 
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Appendix III. Thus the effect of the heat transfer parameters on 
temperature profiles, values of NTU and effectiveness for heat regener- 
ators may be investigated. Results may be compared to the results 
obtained by Creswick [^5 J , Moreland [20j , Bahnke [2j and Nahavandi 
and Weinstein [ 21 ] . However, note the remark at the end of Appendix 
III. 
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APPENDIX I 



PARTICULAR CASE WHERE LONGITUDINAL CONDUCTION IN BOTH THE FLUID AND 
SOLID ARE NEGLECTED. 

The terms describing longitudinal conduction are neglected and the 
differential equations (11) and (12) in the general case are reduced to: 



^ T. 

()ir 

'bXy 



- (>\pY)(u-lL) 

= - 



Let A p V = b 

the Laplace transform of these 2 equations are: 



SIL rr b(\r — U/) 



S _ A(u-ut) 



From equation (29): 



a — 



• tr 



Si- b 



Replacing u. in (30) by its value, we get: 

d\? 



S IT 









or 



lu 






= 0 



(27) 



(28) 



(29) 



(30) 



(31) 
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Solving this ordinary differential equation yields: 



U 



= Gex|o j^-(s+h)y. 



With the boundary condition V (o,s) = J. , we get: 

s 









and 



tt = 



• exb 

S (s + b) ' 






From a table of Laplace transforms, we get the formula 

T 

For a = 1 



L .+ . i,_, [. 1 ^] 



i-e* = I» [a\^] 



( 32 ) 



( 33 ) 



A theorem of Laplace transforms states that if a is any constant and 






= F(s) 



then 



Hence 






( 34 ) 
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From a theorem on the Laplace transform of an integral, we liave 







It then follows that: 




i 

S(S + b) 



o< 



But 



oL 

1 

s(s + b) 




(35) 



(36) 



Hence from (34), (35) and (36), we get: 



r ^ 

= t j bj 



From properties of Bessel functions, we find: 



(37) 



then by integrating by parts, the integral of equation (37) becomes: 











X- 

e 



o 



-h% 






(38) 
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From equations (37) and (38): 









( 39 ) 



The translation theorem of Laplace transforms states that if 

f('t) - H (T:-cu)0(r -a) 

where H — <3u^ is Heaviside's 

unit step function defined by 



H^X-a) - O ^or 'C<C 






H (X — cl) — i ^or 'C' > 



then 







= e 



— C^S 



£j0('c): 



From equation (32), the transform of \r can be written as: 



(40) 



V — 



-H 



Ab 

s 



(41) 



From equations (39), (40) and (41), it is found that: 



0(x) = ^ I ^ = 1 e 

I. O 



40 



and 







then for ^ , we get: 



V — 




I. 









(42) 



Eq (33) can be written as: 



U 



e 




b 

S(s+ b) 



Ab 

e 



(43) 



Using equations (35), (40) and (43), we can write: 



Lt 









(44) 



The equations (42) and (44) were evaluated by the program SCHUMANN. 
Numerical results were checked against solutions from numerical in- 
version of equations (32) and (33) . The agreement was to three 



significant figures . 
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A mor6 simplified solution esn be obtsined by neglecting the term ^ ^ 
in equation (28) which describes the energy stored in the fluid. Thus 
the equations (27) and (28) become: 






Du- 



= - A (v - U/) 



transforming the above equations yields: 



( 45 ) 



(46) 



Sul = b(u-uy) 

— - A ( u- - a) 



(47) 



(48) 



From (47) : 



Lt 



btr 

S + b 



(49) 



Substituting this value of to in (48), we get the ordinary differential 
equation in \T : 



dir 



+ 




\ \j- 

s+ b/ 



The solution is: 

or = C e 



o 
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( 



With 'IT (o,s) = 1 , we find: 



ar — 1 p V 5-*-bl^ 



V = _±_ e 

s 



(50) 



Using equation (39), the inverse transform of (50) can be written as: 



f 

V = 1 1+ 






5 (51) 



From eq. (49) and (50), it follows that: 



iL = e 






u 

£+ b 3- 



S(S + b) 



Using equation (35) yields the inverse transform of (52) 



(52) 



a = €/ 



/ 



- / .-‘■J 



9/^ A b 



cl^ 



(53) 



The solution to this simplified case has been derived by Moreland [ZOj 
in the form of a series. His solution can not be evaluated at y = 0. 

The equation (51) was programmed and checked against Moreland’s solution. 
The agreement is good up to 8 significant figures. 
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APPENDIX II 



PARTICULAR CASE WHERE THE FLUID-SOLID BOUNDARY RESISTANCE IS NEGLIGIBLE, 
I .E. , ha = infinite . 

In this case, ha infinite is equivalent to assuming that the fluid 
and solid temperatures are equal at any time throughout the bed. 

Combining the equations (11) and (12) and letting 'Ll = LL 
yields the following equation: 



(^+-) 
V py y 



'(>\r 



Let c 





( 








(54) 



The transform of equation (54) is: 



CSV = — 

d 

or 

ai^ _ _ csir = 0 

The characteristic equation is then: 

olV^ - r - os = O 
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This is a quadratic equation with complex coefficients, 
be written as : 



r = 




The roots may 



r = 







( 55 ) 



where 






0 

1 

Mac 



Hence : 






( 56 ) 



Applying the boundary condition at y = c>o , we have: 

V ( oO , s) = 0 

Thus only the roots with negative real parts may be used. 
Suppose that r^^ has a negative real part, 

then 
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s. 






Since 



then 



ir (o^s) _ ^ p/'^ 

- s 



2.a ^ 




= e. 


Xe 


s 



From table of Laplace transform we have: 






V' 



er|c 



f 



Since 












then : 
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JL) e 
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X 
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+ e 
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and 






% 



) erfc 


\)^ 










f 



erJ-c 



:t_ 



\[^ 



(57) 



The equation (57) was programmed by using the subroutine ERFN . After 
some testing runs, it was found that for large values of y, the second 
term in equation (57) did not give accurate results for the reason that 
the exponential term becomes very large and the complementary error 
function becomes very small. Thus the product of these two functions 
certainly gives large error. The error can be minimized by approximating 
the complementary error function as a series. A well known asymptotic 
series for erfc x is: 






■f 



ertc X 



( \ L _ Izl. _ + 

V X 5/x* a*- 



) 



Thus the equation (57) may be replaced by: 






v('J^,T)=A^Crfc 









oO 



+ 






(58) 



fit 



m = d 






n-1 2n-l 
OO 
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where 



X = 

(l+ 

^ 

4A(i + )T)((:«y-^1) 

The equation (58) was evaluated by the program JENKINS. Solutions 
from this program were compared with results from the programs TEMFLUl 
and SCHUMANN. Discussions of these results are presented in Section 6. 
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APPENDIX III 



GENERAL CASE APPLIED TO A MODEL OF FINITE LENGTH. 



The following mathematical derivations are applied to heat regener- 
ators and packed beds of finite length. 

The energy balance is as follows: 
a. For the fluid phase; 



Heat stored in an element of fluid : 






Convection by moving fluid : 







Conduction in the fluid : 







Heat transferred to the fluid element 
by convection 



kA (T^-T,) 
L 



then : 





■f + -D X- 






( 59 ) 



(b) For the solid phase : 



Heat gained by an element of solid 









Heat transferred to the solid element by 
convection 



Jlk(Tf -k) 
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1 






I 



Heat transferred by conduction from the 
solid element 




then : 









^ G 









Multiplying (59) and (60) by L : 

ha 







\^A'^ e 





Aka'^ 6 



t<s^s 



_L 7)^Ts. 
hA 






Let us define the following parameters: 



X = 



— 

L 



dimensionless length parameter 



t = 



hA& 

W,c, 



dimensionless time parameter 



A'= kA. 

L 



dimensionless conduction parameter 



NTU = 



hA 



dimensionless heat transfer unit 



Substituting these parameters in (61) and (62) yields: 



(P< 



'DTl 

W.Ci' 'll -b 





-I- 



\ kau'^ X- 




( 60 ) 



Ts) ( 61 ) 

(62) 



) (63) 
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and 



-6 t 



= (^) 






(T,-T.) 



(64) 



Multiplying 


(63) by 




, we get: 


P^A^c^L 


1 

II 


( \l 


( HL 


ki \ 


^ t 




\ hA / ^ 


Vl^ hA/V 



■WsC 

ff 






Let us define: 



oi = JL 

po 



= thermal diffusivity 



P' 









PsC^sAs 

ratio of heat capacities per unit length 



Substituting these parameters in equation (65) yields: 



D h 






If we define: 



It = 



_ Ts 



Ti 



I n4 



XX = 24- 
Ti 
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where is the step input injected fluid temperature, the equations 
(64) and (66) become: 



T> h 



Ct — ^ -h ( U - tt) 

7 ) ^ ^ 



( 67 ) 



"b V 

?) t 



cup> 



/ 7)''ir 



7) X’ 






tc) 



( 68 ) 



where ^ r: 

NTU 

b = -^ 

NTU 

The initial conditions and boundary conditions are assumed as follows: 
(a) Initial conditions : 

The initial fluid and matrix temperature are uniform and equal. 
The base scale temperature can be chosen so that these temperatures can 
be taken as zero for convenience: 

u.(x,o) = \r(x,o) = 0 



(b) Boundary conditions : 

(1) At X = 0 and t = the injected fluid temperature is 
suddenly changed to a different higher value and held constant there- 
after: 



\r (o, -t) = i 
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rises 



(2) At X - 0 and t = 0^, the solid temperature 
to the value of the step input temperature of the 



instantaneously 

fluid; 



= 1 



(3) The matrix is insulated at X = 0: 

■bTs 



^ X 



(o,t) = o 



(4) The matrix is also insulated at X = 1; 

^)T, 



X 



(l,t) = O 



From equation (68) we have: 



tt = 



7 






?)t 






b — + CU IT 



( 69 ) 



7)u 


A. 


“'3V 


as' 4 . 


b 4 




(70) 


dt 














<)^a 




1 

> 

1 




4 


q^7)V 


(71) 


7> x*- 




7)x"^t 


'b x“ 


7i x^ 


?)X"_ 





Substituting equations (69), (70) and (71) in (67) and rearranging the 
terms, we get: 
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0 



(72) 



4 



' 7) X 



/_b \ 




i \ dir 
7 “D-fc 



4 



± 






DV 



The transform of equation (72) is: 




f ab\ d 
' ly / cC 







4 



\ y 4> 




o (73) 



the corresponding auxiliary equation is then; 



A,r 



4 A^r = o 



(74) 



Where > > complex coefficients of equation (73). 

The general solution in the Laplace S plane for the fluid temperature is: 



V = CJs)e/'\ C,(s)e^‘+ C3(s)e"‘V 



(75) 



where r]^, T 2 , ^3> ^4 complex roots of equation (74). The boundary 

conditions are transformed and then used to determine the coefficients 

BC.l V (o,s) = i 

s 

BC . 2 iL (o, s) = ^ 

s 
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BC.3 


-^(o.s) . 0 
*l)x 




BC.4 


(l.s) = 0 




Applying BC. 1 


to equation (75) gives: 




v(o,s) 


= rc„(s) = j_ 

^ ^ s 


(76) 



Taking the derivatives of XT (x,s) with respect to x and evaluating 
them at X = 0 yields: 



l)\r 

X 



(o.s') 



YZ 

m = i 



(77) 



()’-U 



(o,s) 



i 



T) 




Y}n c« 



(78) 



(79) 



From equation (69) we have: 



Z)x L 'bx ^ ?)x» ^x^ ^'axj 



(80) 



Applying the BC (3) to equation (80) and using the equations (77), (78) 
and (79) give: 



k k h ^ 

mszl. 



( 81 ) 



55 



Applying the BC 4 resulting in: 



- 






7}zz t 



i 'V%mi 



r . C , 



Let us define the quantity: 



Rn = 




')^n + brj + (9+S ) r^ 













then the equations (81) and (82) can be written as follows: 



4 

D 



= 0 



and 




= 0 



Finally, applying the BC (2) to the transform of u yields: 



u(o,s) = i 



_a ^ ^ (o,s) H- b 

I ^ 



or 



or 



EVn = - 



i 



where = 



•a A'' r^ + br„ 
p n n 



•Dir 

^x 

JL 



+(s+ ij; )v(o,s) 



= 

-Tl » 4. 



= -1 



= o ( 82 ) 



(81a) 



(82a) 



1 (83) 



(83a) 
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The equations ( 76 ) , ( 81 a) , ( 82 a) and ( 83 a) are then used to solve for 
the coefficients C^(s) of equation (75). 

We have the following matrix equation; 



1111 




Cj(s) 




1 

s 


Rl R2 R3 R4 




€2(5) 




0 


R^e^^ R26^^ R3e^3 R^e^^ 




C3 ( s) 




0 


Zi Z2 Z3 Z4 




C4(s) 




-1 



( 84 ) 



the determinant of the matrix ( 84 ) can be written as: 

/\ = R 2 R 3 (Z 4 -Zi)(e’^ 3 _gr 2 ) + R 2 R 4 (Z 3 - Zi)(e’^2 . gr 4 ) 

+ R3R4(Z2-Zi)(e’^^-e’^^) + RiR3(Z4 - Z2)(erl - 

+ R3^R4(Z2-Z3>(e’^^-e’^^) + R^R^(Z4 - Z2)(e^^ - 



the coefficients are then: 



Ci(s) = 1 



= 1 
A 



1 1 
s 



r 2 „ r 3 



0 R2C R 3 e R4e 



,r 4 



-1 



1 

s 



|^R2R3Z4 (e’'^ . + R2R4Z3(e'"2 . + R3R4Z2(ei^^ - e’^3)j 

j^-R3R4(e’'^ - e’'^) + R2R4(e’"^ - e’^^) _ R^R3(e’'3 _ e’"^)j j 
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tf 




Co(s) = 1 

A 



1 

RiG 



rl 



I 

s 

0 

0 

-1 



1 

R3 

Roe 



I 

>^3 



i \ -i 

A 1 s 



RlR3Z4(e^3 . ^ RiR^Z3(e’^l - + R3R4Zi(e’^'* - e>^ 3 )J 



.R3R4(e>^^ - e’^3) . R^R^(erl . ^ R^R3(erl . grS) 



( 86 ) 



Co(s) = 1 . 

A 



1 

Rl 



rl 



R2e’^2 0 



1 

R^e 



r4 



-1 



^ I i |^RiR2Z4(e’^2 . ^rl) + R^R4Z2(e’^l - e’^^)+R2R4Zi - e^^)j 
+ jjR2f^4(f'^^ ■ + R]^R4(e’^^ - + R]^R2(e’^2 . | 



(87) 



C4(s) = 



1 

£ • 



1 

^1 

Rie 



rl 



1 1 

^2 ^3 

R2e"^ R3e’^3 



_1 

s 

0 

0 

-1 



^ I - 1 jR]^R2Z3(e’^^ - e’^^) + RiR3Z2(e’^^ - e’^3) + R2R3Zi(e’^3 . 



j^-R2R3(e’^3 _ gr2) + Rj^R3(e’^3 _ ^rl) _ r^R2(c 






( 88 ) 
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r> 




A subroutine may be written for the equations (69), (73), (75), 
(85-88), Numerical results for temperature profiles may be obtained 
by using this subroutine with the main program of TEMFLUl. If there 
are difficulties with the solution to this case, such difficulties 
probably have their source in the assumed boundary conditions. It 
seems that the boundary conditions (2) and (3) are not independent. 
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APPENDIX IV 



PROGRAM LISTINGS 



1 . PROGRAM TEMFLUl 

a. PURPOSE : 

This program finds the inverse transform of v and u, using 
Salzer*s method of numerical inversion. 

b. USAGE : 

(1) INPUT FORMATS 

The input data are read from two cards. The first card 
reads 8 parameters in floating point format 8F10.5. The second card 
reads 3 parameters in floating point format 3F10.5 and the run number 
M is fixed point format 13: 

TKS = Solid thermal conductivity, 

TKW = Fluid thermal conductivity, 

Density of solid phase. 

Density of fluid phase. 

Specific heat of solid phase. 

Specific heat of fluid phase. 

Porosity of porous media. 



ROS 

ROW 

CS 

CW 

POR 

HA 



Btu/hr 

tt 

lb mass/ft^ 
lb mass/ft^ 
Btu/lb mass ^ 

ft 

dimensionless 



Heat transfer coefficient, based on a 
unit volume of bulk porous media. 



VEL = Fluid interstitial velocity, 

XI = Distance from point of fluid injection, 

F = Number of time units per hour, 

M = Run number - Set M=0 on last data card 
to stop the program 



Btu/hr ff^ OF 

ft/hr 

ft 

time units/hr 
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1 




A 

B 

C 

Y 
T 
N 

V 
U 

ERV 

ERU 



c , 



(2) OUTPUT FORMATS 

= Ratio of thermal dif f usivi tics , 

= Ratio of thermal conductivities, 

= Dimensionless parameter 

= Dimensionless distance 

= Dimensionless time 

= Order of polynomial 

= Fluid temperature fraction 

= Solid temperature fraction 

= Difference between two values of V using 

two adjacent orders of polynomial 

= Difference between two values of U using two 

orders of polynomial 

SPECIAL INSTRUCTIONS 



The program TEMFLUl calls the subroutine VUBARl 
the function AITKENF to interpolate V and U. 



d. MATHEMATICAL METHOD 

See section 3 and 4 of this thesis. 



dimensionless 

It 



adjacent 



It also uses 
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PROGRAM TEMLUl 



START 
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( 



r 




G=XR(N,J)/T 

W=XI(N,J)/T 




TEMPV=TEMPV+TVR 

TEMPU=TEMPU+TUR 




VB(N)=TEMPV/T 

UB(N)=TEMPU/T 



COMPUTE ERROR 

ERV(N) ®VBCn)-V&(n-1) 
= UB(n)-UBCn-I) 
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): 





DO 700 K=1,I 



FV(K)=TV(K) 

FU(K)=TU(K) 




V=AITKENF(Z1,FV,X,I-1) 
U=AITKENF (Z 1 , FU , X , I - 1 ) 




© 



1 
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I 



ll 



2 . 



SUBROUTINE VUBARl 



a. PURPOSE : 

This subroutine, given the values of the dimensionless para- 
meters and of the transformed variable S, calculates the complex co- 
efficients of the quartic equation (18) , calls the subroutines POLYRT 
or COMSUB that find the complex roots of equation (18), calls the sub- 
routine POLYVAL to check the accuracy of the roots, then selects the 
roots with negative real parts to calculate the coefficients and C 2 
of equation (19a) and finally computes VR, VI, UR, UI. 

b. USAGE : 

(1) INPUT ARGUMENTS ; 

G = Real part of the transformed variable S 
W = Imaginary part of the transformed variables 
A = Ratio of thermal dif fusivities 
B = Ratio of thermal conductivities 
C = Dimensionless parameter 
Y = Dimensionless distance 

(2) OUTPUT FORMATS : 

(a) **N = ,J = ,MG = , PZR OR PZI IS LARGER THAN l.E-4** 

is printed if the roots are not accurate. ”N** is the order of polynomial; 
”J** identifies one of the zeros of this order; **PZR** and **PZI** are the 
values of the quartic equation (18) evaluated at the root. **MG** refers 
to the subroutine used for solving the quartic equation; **MG = 0** or 
'*MG = 3** is printed if POLYRT has been used; "MG = 1** or "MG = 2" is 
printed if COMSUB has been used; "MG = 4" or "MG = 5" is printed if both 
subroutines have been used; the number is 4 if POLYRT has been used first, 
and 5 if COMSUB has been used first. 
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(b) "N = ,J = ,MG = , ONE ROOT HAS NEGATIVE REAL PART** 

is printed if only one root with negative real part has been found. 

(c) '*N = ,J = ,MG = , THREE ROOTS HAVE NEGATIVE REAL 

PART” is printed if three roots with negative real part have been found. 

(d) ”N = ,J = ,MG = , CONSTANT VECTOR NOT EQUAL TO CHECK 

VECTOR” is printed if the calculation of C^ and C2 has not been accurate. 

(e) VR = Real part of the transform of v 

(f) VI = Imaginary part of the transform of v 

(g) UR = Real part of the transform of u 

(h) UI = Imaginary part of the transform of u 

(i) VR, VI, UR and VI are set equal to zero if one of the 
outputs (1), (2), (3) or (4) is printed. 

c. SPECIAL INSTRUCTIONS : 

VUBARl uses the subroutine MULT for multiplication of two 
complex numbers. 
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SUBROUTINE VUBARl 
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I 

i 
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137 
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DO 200 1=1, A 
RR(I)=0 
200 RI(I)=0 




PRINT 

PZR OR PZI IS 
LARGER THAN l.E-A 



GO TO 110 
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3. 



PROGRAM JENKINS 



a. PURPOSE : 

This program finds the solution to the special case where the 
fluid and solid temperature are assumed to be equal. The equation (58) 
is programmed, using the function ERFN which calculates the error 
function. 



INPUT FORMATS: 



A 


= 


Ratio of thermal dif fusivi ties , 


dimensionless 


B 


= 


Ratio of thermal conductivities, 


M 


C 


= 


Dimensionless parameter , 


M 


T 


= 


Dimensionless time 


It 


M 


= 


Run number. Set M = 0 on last data card to 
stop the program. 




c . OUTPUT FORMATS 




X 


= 


Dimensionless distance 




ERC 


= 


Value of the first term of equation (58) 




E2 


= 


Value of the second term of equation (58) 




V 


= 


Fluid temperature fraction 




and E2 


are 


printed out to show their relative importance 


• 
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PROGRAM JENKINS 



100 READ 
A,B,C,T,M 







75 




COMPUTE SERIES 
OF ERFC(y2) 



F=l. 

SUM=0 
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4. 



PROGRAM SCHUMANN 



a. PURPOSE : 

This program finds the solution to the special case where both 
longitudinal conduction in fluid phase and solid phase are neglected. 

The equations (42) and (44) are programmed, using the subroutine 
GAUSSN to evaluate the integrals. GAUSSN itself calls the subroutine 
FOFX which evaluates the integrands by using the subroutine BESSELL to 
find the values of modified BESSEL functions of first kind (order 0 and 1) . 

b. USAGE : 

' The input and output format are the same as those defined in 

program Jenkins. 
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PROGRAM SCHUMANN 



1N1T=-1 




Y=Y+10 . 
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